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Abstract. We explore a series expansion method to calculate 
the modes of oscillations for a variety of uniformly rotating 
finite disks, either with or without a dark halo. Since all mod- 
els have the same potential, this survey focuses on the role of 
the distribution function in stability analyses. We show that 
the stability behaviour is greatly influenced by the structure of 
the unperturbed distribution, particularly by its energy depen- 
dence. In addition we find that uniformly rotating disks with 
a halo in general can feature spiral-like instabilities. 
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The beautiful spiral structure, bars and rings which disk 
galaxies often show point at the presence of (possibly tran- 
sient) perturbations in the disk. The visual appearance proba- 
' bly should be explained to some extent by gas dynamics, but 
the stellar component must be an important player too when 
\Q the spiral structure is well-developed and thus probably not 
too ephemeral. 

The question of the stability of such disks and the structure 
in of the instabilities which might occur is an important one, but 
5 the solution turns out to be quite involved for realistic cases. 
*~^~} Currently, N-body calculations offer the most flexible tool for 
the study of the evolution of models for such galaxies. Alter- 
C natively, one models the perturbations in a galaxy by search- 
es [ing for linear modes, i.e. selfconsistent linear solutions of both 



the Poisson equation and the linearised collisionless Boltzmann 



I 

2 equation. The results from this approach seem to be compatible 
+^ with N-body simulations, at least concerning the early stage of 

_j the instability (see e.g. Sellwood & Athanassoula, 1986). Al- 
though the literature offers general methods for calculating 
these modes (Kalnajs, 1977), up to now the calculations only 
have been done for a few isolated cases. 

In this article, we calculate linear modes for a particu- 
lar set of simple models, i.e. flat disks with finite radius in 
a quadratic potential. Although these models are an oversim- 
plification since they do not show differential rotation, they 
can reveal some interesting information about the stability of 
real galaxies, particularly in the central parts where the po- 
tential is sometimes approximately quadratic. We have chosen 
these models since we can apply a general and straightfor- 
ward scheme to calculate the linear modes in quadratic po- 
tentials. This method is based on the series expansion approx- 



imation for the solution of the linearized Boltzmann equation 
(Vauterin & Dejonghe, 1995), which becomes exact in the case 
of quadratic potentials. 

The mode analysis has already been performed analytically 
for uniformly rotating self-consistent disks without a halo and 
with some particular distributions (Kalnajs, 1972). We extend 
this work by calculating modes for a variety of different disk 
distribution functions, embedded in an inert halo. 

1. Method 

1.1. Solutions for the linearized Boltzmann equation 

The total disk is supposed to consist of two parts, a time- 
independent axisymmetric unperturbed part, plus a small per- 
turbation. The distribution function is denoted as 



f(r, 8, v r , ve, i) = f (r, v r , ve) + f'(r, 8, v r , ve, t), 



(1) 



while the total potential, which is defined as a binding energy 
(decreasing outwards), reads 



V(r,8,t) = V (r) + V'(r,8,t). 



(2) 



The unperturbed quantities are marked with a subscript 0. 
Moreover, using Jeans' theorem, we know that the unperturbed 
distribution function can also be written in terms of the two 
integrals of motion, the binding energy E and the angular mo- 
mentum J, defined by 



E = V {r)- 1 -{vl+v 2 e ) 



and 



J : 



rvg. 



(3) 



(4) 



Using Poisson brackets, the linearized collisionless Boltz- 
mann equation for small perturbations is quite compact: 



[fo,V] 



(5) 



This fundamental equation describes how the motion of the 
stars in a given perturbing potential determines the evolution 
of the distribution function, up to the first order in the pertur- 
bation. 

The right hand side of the equation can be written as 
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(6) 



pendence fo(E, J), this becomes 

" (§§ (jB ' j)WvE + ¥f^ E ' j)v ^ j ) ■ VrV '' 

which can be transformed again into 
^[E,V'] + ^[J,V]. 



dE' 



d.J 



(7) 



(8) 



The operator acting on /' in the left hand side of the equation 
(5) is linear and is "seeing" E and J as constants since these 
are integrals of the motion. Therefore we can write the solution 
as 



ri _ dfo dfo 
1 ~ dE lE dJ Jj ' 



(9) 



with f' E and f' } the solutions of the linearized Boltzmann equa- 
tion having respectively fo(E, J) = E and fo(E, J) = J as 
unperturbed distribution function. For a given unperturbed 
potential V(r), one only has to solve the Boltzmann equation 
in these two special cases and use the linear combination to 
calculate more general distribution functions. 

In this paper, we assume a quadratic potential: 



V (r) 



1 2 2 

- — il r , 

2 ' 



(10) 



with fio the vibration frequency of the stars in rectangular 
coordinates. We will use the power series expansion method 
(Vauterin & Dejonghe, 1995, hereafter paper I) to calculate 
the perturbed distribution. In this approach, the radial coor- 
dinate and the velocities are expanded in power series, while 
the angular dependence and the time is expanded harmonically 
(the circular velocities are chosen as zero points for the veloc- 
ity expansion). In the linearized case, the different harmonics 
are independent and so the expansions can be written down 
for each individual harmonic m. Adopting the same notations 
as in paper I, this reads 



/' 
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-mff) 
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and 



T rt \ ^ m i(u:t- 

V = y ^a k e y 

k>0 



mff) 



(11) 



(12) 



(13) 



with co/m the pattern speed of the perturbation. When these 
forms are substituted in the linearized Boltzmann equation, 
one can collect terms with equal powers and construct a set 
of equations in order to determine the unknown vT,j,k ( see sec ~ 
tion 2.2 of paper I). This set of equations is represented most 
conveniently in matrix notation, by defining the vectors 



(pT.n,, 



( Po,n,k ) 



Pl,n — l,k ) 



Pl,n- 



Pn — l,n,k ) Pn.O.k ) : 



Pn 



Pn,0,k 



(14) 
(15) 



In addition, we define three sets of matrices (they are defined 
by enumerating the nonzero elements) 



[A n (v, s)]( n+ l) x (n+l) 



Ak,k+i = sk 

A n—k+l 

Ak + i,k = 



(16) 



+ 1) Xn 



Bk,k 


= m 




Bk,k + 1 


= ki 


(17) 


Bk + l,k 


= i(p — n + k) 





and 



n+l) x(n-\-2) 



Dk,k 
Dk,k + 1 



m(n + 2 — k) 
ikp 



(18) 



If we further rescale the perturbation pattern speed into a new 
parameter 



2Q 



(19) 



the resulting set of equations for the unknown P™ p reads in 
the case of a quadratic potential 



2fio*4n(-fo, i)P™ p 

*2 m ~p m \ ^ m-rxm pO 

*^n,p+l — l,p+l / J &q "n,q*n-\- 1 ,p — q-\- 1 • 



(20) 



q=0 



The matrix A„( — vo,i) is regular for all non-integer values of 
vo- The integer values correspond to resonances. In the non- 
resonant case, this equation can easily be solved recursively. 
When E and J are used as unperturbed distributions (using 
(9), the response distributions of these simple functions are 
sufficient to calculate a general solution), the expansion of the 
unperturbed distribution contains only a few low order terms. 
It is then sufficient to choose a perturbing potential having a 
polynomial form in order to obtain a terminating expansion of 
the perturbed distribution. In this case, the method becomes 
exact. The total degree of the response is the sum of the maxi- 
mum degree of the unperturbed distribution and the degree of 
the polynomial perturbing potential. 

1.2. Solutions for the Poisson equation 

Since we suppose that the unperturbed disk extends only up 
to a maximum radius d, it is reasonable to use the family of 
potential-density pairs described by Hunter (1963): 



PHun.i 



1 9 r p ™ W^w) , ut _ 
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(21) 
(22) 



with m > 0, I > |m|, Pi, m the associated Legendre functions 
and 



9i 



r(^ + i)r(^ + i) 
r ( 1± f ±i ) r ( i =f ±i ) 



(23) 



For finite disks, these densities present a complete set and can 
be used as a basis for the expansion of general perturbed den- 
sities. 



The construction of self-consistent modes for a particular dis- 
tribution and with a particular harmonic symmetry m consists 
in the search for solutions /' of both the Boltzmann and the 
Poisson equations. 

We write a general perturbing potential as a linear combi- 
nation of (22), 



V 



n ,m-\-2p 



(r,e,t). 



(24) 



p=0 



In principle, s should be +00. In practice, one has to set a rea- 
sonable upper limit for s without compromising on accuracy. 
This can easily be checked afterwards. Typically, s ~ 12 in our 
calculations. 

Since the potentials (22) are polynomials in r, the series 
expansion method can be used to calculate the response per- 
turbed distribution function for each of them. The result will 
of course depend on the pattern speed lo. In addition, the in- 
tegration over the velocities yields the perturbed mass density 



-trm r m / \ m / \ 

v mm,m + 2p Jm + 2p\ ul ) Pm + 2p\ u> ')■ 



(25) 



Each of the perturbed densities, which are also polynomials, 
can be expanded in the densities (21), again up to a maximum 
order s 



Pm + 2p 



Cq,p(^)PHu 



n ,m-\-2q 



(26) 



For our models, these coefficients can be calculated analytically 
using the expansion formulae given by (Hunter, 1963). 

Using the matrix c qiP , one can now write the response po- 
tential resulting from the perturbation (24) as 



X 



p=0 



c q ,p(ui)a p 



Hun ,m-\-2q • 



(27) 



For a self-consistent mode, this response is required to be equal 
to the original perturbing potential. This can only be achieved 
for particular values of the pattern rotation speed lo. For gen- 
eral values of lo, one can require the response potential only to 
be proportional to the perturbing one, with a scale factor A. 
This results in an eigenvalue problem for the (s + 1) x (s + 1) 
response matrix C(lo), with elements c q , p (w), yielding eigen- 
vectors a p and eigenvalues A which are the above mentioned 
scale factors. Note that in the present case, the elements of the 
matrix C are found exactly and in a form that keeps lo analytic. 
However, since C has an infinite dimension and is not diagonal, 
the eigenvalues can be calculated only approximately. 

1.4- The eigenvalue problem for uniformly rotating potentials 

One finally has to calculate the eigenvalues for all values of lo in 
the complex plane (we will call this relation the eigenvalue re- 
lation A(w), which is multiple- valued). Values for lo resulting in 
a unity eigenvalue correspond to normal modes. Generally, the 
search for unity eigenvalues has to be performed numerically. 
Normal modes with a complex lo occur in complex conjugate 
pairs, and one of them will have an exponentially growing am- 
plitude, causing an unstable behaviour. On the other hand, real 



without any exponential growth. We call these modes stable. 

In a quadratic potential, the response behaviour is partic- 
ulary simple since all resonances are global: when the pertur- 
bation rotates with a pattern speed to such that lo/CIo has an 
integer value with the same parity as m, the whole disk is in 
resonance, while otherwise no part of the stellar distribution is 
resonating. This is a direct consequence of the fact that, in a 
harmonic potential, all stars have exactly the same vibration 
frequency 2fin. So the eigenvalue relation X(lo) has poles at 
integer values of lo/CIo with the same parity as m (see fig. (1) 
for an example). 
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Fig. 1. Eigenvalue relation X(to) along the real axis for model I with 
Q = 0.5, a = 1 and m = 2. Each curve corresponds to a branch of 
the characteristic equation that yields the eigenvalues. 

In addition, the functional dependence of the elements of 
the response matrix c qiP (Lo) is particulary simple, since they 
consist only of a sum of fixed poles: 



r,pH = X/ 



Qol 



(28) 



with a summation index I having the same parity as m. This 
is easily shown in the framework of fourier expansions along 
the orbit (e.g. Kalnajs, 1977), since all orbits have one single 
radial period 2fin. 

The search for unity eigenvalues (corresponding to a normal 
mode) is equivalent to solving the equation 



\C(lo) 



I\ 



0, 



(29) 



where I is the (s + 1) x (s + 1) unity matrix and |...| stands 
for the determinant. When this equation is multiplied enough 
times with factors lo — Qol to remove all poles, it becomes a 
polynomial function in to with real coefficients (note that these 
factors can never become zero for physical responses since this 
would correspond to singular resonances). The degree of this 
polynomial is constant for a certain value of s (apart from acci- 
dental cancelation of some terms, which exceptionally occurs). 
For distributions which have symmetrical J dependence with 
respect to J = 0, this degree equals s(s + l)/2 + sm, while for 
other ones the polynomial has s more roots. 

The modal frequencies can now easily be found as the roots 
of a polynomial. This is a big computational advantage, as it is 
not longer necessary to perform a brute-force numerical search 



from this, it is sufficient to observe X(ui) along the real axis to 
have a complete view of the stability behaviour. The degree 
of the polynomial gives the number of modes. Every pair of 
modes which is "missing" on the real axis should occur in the 
complex plane, with opposite signs for the imaginary values. 

Between two poles, there are three possible situations con- 
cerning X(ui) (this relation is actually multiple-valued, but 
holds for each branch). 

— The X(ui) approaches the asymptotes with opposite sign. In 
this case, there will always be a real mode inbetween both 
asymptotes. 

— The X(ui) approaches both asymptotes with negative sign. 
There can be a pair of real modes if X(ui) reaches a value 
greater that one. We never observed this though: in all our 
models A(w) remained negative between such poles. 1 

— The X(ui) approaches both asymptotes with positive sign. 
In this case the curve X(ui) may or may not intersect the 
A = I line. The corresponding frequencies are real resp. 
complex, and the mode is stable resp. unstable, depending 
on the details of the unperturbed distribution function. 

In the first two cases, all modes are stable. In the last case 
however, unstabilities can occur for particular unperturbed dis- 
tributions. This is why we will call this behaviour "potentially 
unstable" . 

2. The models 

In the following paragraphs, we will model a galaxy using two 
components: 

- A disk, which can be subject to perturbations. This disk 
has a central mass density Io and extends up to a radius d. 
The distribution function is written as fo(E, J), while the mass 
density is represented by po(r)8(z) (with 8 the Dirac function). 
This disk creates a potential which is denoted by Vr>,o(r). In 
general, this potential is not quadratic. 

- A spherical halo, which is supposed to be unresponsive 
to perturbations in the disk. We represent the mass density of 
this halo by pjj(r), while the potential is denoted as Vh,q(t). 

The total potential is quadratic 



H/D 
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V (r) = VD, (r) + Vir,o(r) 



1 n ^ 



(30) 



In practice, we choose a particular /o, defining a disk mass 
density po(r) by integration over all velocities. This mass den- 
sity produces a potential Vd,o which is calculated using the set 
of potential-density pairs (21, 22) (only m = is used). Then 
the spherical halo density pu is adjusted in order to produce 
an overall potential Vo of the form (30): 



Ph{t) 
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r 2 dr 



: dV Dil 
dr 



(31) 



which expresses the spherical halo mass density for all values 
of r between and d. Further away the halo mass density can 
be chosen freely. Of course, pu should not be negative for any 
value of r between and d. For a given Vr>,o and fio, this yields 
an upper limit for the disk mass density. The halo strength will 
be quantified by calculating the ratio of the total mass in the 
halo in a sphere with radius d to the total disk mass: 



J It follows from the sequel that such modes would cause a disk 
to become unstable with increasing halo mass. 



M d 2tt J Q d p (r)rdr 



(32) 



2.1. Rotating isotropic models with decreasing f(E* t ) 

The algorithm proposed in the previous paragraph calculates 
the self-consistent modes for the disk having a general unper- 
turbed distribution function fn,o(E, J). Since we want to dis- 
cuss how the stability of the models depends on a few impor- 
tant parameters of the distribution function, we will start with 
a further simplification by focusing on functions which are of 
a form similar to those used by (Kalnajs, 1972). In order to 
simplify the notations, we define 



E* n = i(Qo - Ut 2 )d 2 + E + QJ, 



(33) 



with fi a constant parameter which lies between — fio and fio. 
This can be written as 

E* n = - fi 2 )(<2 2 - r 2 ) - i[(„ e - fir) 2 + „*]. (34) 



2 



2 L 



In the following sections, we will set the unit of distance equal 
to d, so that the disk extends up to radius I. We will now use 
distribution functions of the form 



fD,o(E, J) 



/w E* n > 
EX < 0, 



(35) 



with real powers ct t > —1/2. These models are isotropic in a 
frame rotating with a speed fi (note that outer edge of the dis- 
tribution is not given by the escape velocity, but by E* t = 0, 
which is a circle centered at v r = 0, vg = fir). For this reason 
we call them "rotating isotropic". The distribution has a dis- 
continuity at its edge if there occur terms with ct t < 0. Special 
attention should be paid to distributions which have a discon- 
tinuity in their functional dependence, since this gives a Dirac 
^-function singularity in the response distribution. This singu- 
larity is integrable though, and adds only a finite contribution 
to the response mass density. 

The Kalnajs disks are self-consistent models with a 
quadratic potential and with a distribution which is isotropic 
in a rotating frame. They are one of the standard models con- 
cerning stability of disks in literature since all results are ana- 
lytical (Kalnajs, 1972). However, they have an integrable, but 
very unphysical, singularity at the edges. In addition, they have 
an "inverted" distribution function, i.e. the number of stars is 
decreasing for increasing binding energy. In a real galaxy, this 
will usually be the other way round. If the stability analysis 
of the spherical case is of any guidance here, one might expect 
that the slope of the energy dependence is an important factor 
for stability (Fridman & Polyachenko, 1984). 

For this reason, we analysed some models for which the 
number of stars is everywhere a non-decreasing function of the 
binding energy. We have chosen a very simple set of models, 
having 

dfo,o = f1i[E* n ] a , (36) 
with a positive and integer. The mass density is given by 



(37) 



(Q2_ Q2)(l _ r2); 



a + 2 V 2 
and the streaming velocity reads as 



(ve) 



Or. 



(38) 



(39) 



This velocity dispersion is proportional to fig — For a 
constant On, increasing Q means increasing streaming velocity 
and decreasing velocity dispersion. For Q = On, the disk has no 
random motion and is supported by rotation only. In the next 
sections, we will call the distributions given by (36) "Model I". 
These disks usually require a substantial amount of halo mass 
in order to produce a quadratic potential. Even for the case 
closest to self-consistency, a = 0, H / D has a minimum value 
of 2.25, which means that the inert spherical halo is more than 
twice as heavy as the disk. 




Fig. 2. Disk (full lines, as surface density) and halo (dashed lines, 
as volume density) mass density with a = 1 for different values of 
H/D. The curves with highest disk central masses correspond to the 
curves with lowest halo central masses. If no disk is present, the halo 
has a constant density po = 3Qq / 4ttG. 

In addition, we can adjust the halo/disk proportion H/D. 
Fig. 2 shows the mass density for disk and halo for different 
values of H/D in the a = 1 case. Note that the amount of inert 
material increases in the outer regions. The lowest value which 
is consistent with an everywhere non-negative halo density is 
H/D = 6.645. This value is so high since we have chosen a 
spherical halo. For an oblate or even a disky halo, this minimal 
value is much lower. The common intersection point for all halo 
mass density curves corresponds to the point where 



2 dVn.o , x 
r — ; (r) 



dr 



(40) 



reaches a maximum. 

The real axis cut of A(w) for a = 1, Q = 0.5 and m = 2 
is shown in Fig. 1 in a minimal halo configuration. In this 
case, there are no unstable modes, but there exist pairs of 
potentially unstable modes between lo = and lo = 2. The 
qualitative behaviour of X(lo) is quite representative for other 
cases of Model I with different a, Q or m. There is always one 
single region between two adjacent poles where all models have 
pairs of potentially unstable modes. These potentially unstable 



m = they occur at lo = 0), in other words they are corotating 
with the streaming velocity of the galaxy. 

The minimum reached by the eigenvalues in this potentially 
unstable region is increasing for increasing Q and for decreas- 
ing halo-to-disk fraction H/D. If this minimum is greater than 
1, the pair of real solutions will become a pair of complex con- 
jugate solutions, corresponding with instabilities. For a = 
and m = 2, fig. 3 shows the positions of the stable and un- 
stable modes for varying Q in the minimal halo configuration. 
Although the a = case is a somewhat unphysical limiting 
case, we will present most of the results for this model since 
the instabilities appear at the lowest values for Q, resulting in 
clearer graphs. However, it is important to note that models 
with higher a show the same behaviour. Even combinations of 
terms with different a show similar behaviour as long as the 
distribution is an increasing function of the binding energy. 
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Fig. 3. Mode frequency diagram for Model I with a = for m = 2 

2.1.1. Stable modes 

For the a = 1 configuration, Fig. 4 shows a few density profiles 
of the lowest order oscillating (stable) modes, resulting from 
the normal mode calculations. For each value of the symmetry 
number m, multiple modes are found. These modes are further 
catalogued using a second parameter n, starting at 0. For in- 
creasing n, the number of zero-points in the interval [0, d] is 
increasing. Each mode (m, n) has zero points interleaving the 



zero points does not occur. This is reasonable since the total 
mass contained in each mode should vanish (the disk cannot 
lose or gain mass). For m = 0, it follows immediately that the 
radial dependence should have at least one zero point. In gen- 
eral, each mode (m, n) appears to occur at 2n + 2 values for to. 
For small values of r, every mode (m, n) behaves like r m (see 
also paper I). 

Only the lowest order modes have an easy physical inter- 
pretation. The (0,1) mode is a so-called "breathing mode", 
indicating that the galaxy is periodically expanding and com- 
pressing. In the (1,0) mode, the centre of the galaxy is dis- 
placed and rotates around the centre of the unperturbed po- 
tential. The (2, 0) mode essentially produces a rotating bar. 
When no halo is present, the (1,0) mode does not rotate and 
represents a simple displacement of the system. 



m= o 




m= 2 



0.0 Radius 1.0 

Fig. 4. Density profiles of the modes with lowest orders for Model 
I with a = 1 



2.1.2. Unstable modes 

When Q is large enough and H / D is low enough, i.e. when the 
disk is cool and heavy, some of the oscillating modes become 
unstable. Fig. 5 shows the instability limits for the lowest or- 
der modes as a function of the rotation parameter Q and the 
halo-to-disk proportion. The lower right corner contains the 
unstable configurations. This corresponds to cold, rapidly ro- 
tating disks and light halo's. The first instability appears at 
fi/fio = 0.81. At this point, the value of Toomre's stability 
criterion in the centre of the disk is Q = 1.88, lying quite close 
to the value 1.68 in the case of the uniformly rotating sheet 
(Toomre, 1964). 

Fig. 5 further shows that the stability increases rapidly with 
increasing n. For n > 3, there are no unstable modes found. 
However, for unstable modes, the number n loses some of its 
meaning since it does not count the number of zero points in 
the radial density profile anymore. In fact, for values of Q very 
close to fin, even an n = mode contains multiple zero points, 
resulting in a somewhat patchy perturbation. The second pa- 
rameter n is therefore only defined as a continuation of the 
stable case. 

From fig. 5, it is clear that the slopes of the instability lines 
are increasing as m increases. As a result, for large values of 
fi (meaning cooler disks), the disk becomes more unstable for 
larger m. This feature, combined with the remarks in the pre- 



(hot disks), low order global modes (m = 0, 1, 2) dominate the 
instability space. Contrarily, for high values of fi (cool disks), 
short ranged local instabilities dominate (cfr. the Jeans insta- 
bility). 

As can be seen from fig. 3, the pattern speed of the in- 
stabilities increases as they become more and more unstable 
(higher fi), but it remains smaller than the rotation speed of 
the stars. 



m = 1 
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Fig. 5. Instability limits for Model I, a = (lower right corner= 
unstable region). Note that the n = 0, m = modes do not occur, 
as explained in the text. For n = 2, the m = mode is already 
completely stable for all Q. 

Fig. 6 shows the (2,0) mode, which is the most unstable one 
for intermediate values of fi and shows a spiral-like behaviour. 
Most of the unstable modes show a spiral-like structure and 
they are always trailing. 




Fig. 6. Mass density of the unstable (2,0) model I for a = and 
Q/Qq = 0.85. The streaming velocity goes counterclockwise . 

The presence of spiral modes is usually associated with dif- 
ferential rotation and the fact that the perturbation is winded 
in the neighbourhood of resonances. It is clear that this effect 



these have no differential rotation. However, the distribution 
function itself can do the job since it is able to put, depend- 
ing on the radius, varying weights on the various resonances. 
Suppose that the disk is perturbed by a rotating quadrupole 
potential which is aligned to the X-axis of the rotating frame: 



V' = Re [V'(r)e 



i(20-(uv + u, c )t)l 



(41) 



with V'(r) a real function. The perturbed mass density is given 
as a sum of contributions from different resonances, each hav- 
ing the form of a pole (see also eq. 28) 



p > = Re f Al 1 

lUO r + IU0 C — 2ffin J 



(42) 



If Lo c is not zero, the contribution coming from such a resonance 
is not aligned to the X-axis anymore: 
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-[(cor - 2IQ )cos26 l + w c sin20], (43) 



resulting in a response which lags behind the perturbing po- 
tential by an angle a given by 



tan 2a 



2lfin 



(44) 



From fig. 3, it is clear that all instabilities occur slightly below 
corotation (for the whole disk). This means that the lag angle 
is large (and positive) for CR and much lower for ILR, OLR 
and higher order resonances (as 2lfin — u r is smallest for CR). 

In addition, the velocity dispersion is decreasing with in- 
creasing radius, as in most cases. This means that the outer re- 
gions are biased in favour of circular orbits, supporting mainly 
corotation resonance, while the central part has more eccentric 
orbits, giving more contribution to higher order resonances. 
Therefore, it can be inferred that the lag angle is increasing as 
r increases, resulting in a trailing spiral structure. 

Of course, this scheme does not conflict with the existing 
explanation for the formation of trailing spirals, but it rather 
presents an extra effect which is presumably only dominant in 
regions with almost no differential rotation. 

2.2. Rotating isotropic models with increasing f(E* l ) 

As mentioned before, the mode analysis for the Kalnajs 
disks was already performed fully analytically earlier 
(Kalnajs, 1972). These disk have an unperturbed distribution 
function 



fo = E* t 



-1/2 



(45) 



In this special case (and for all weighted integrals of (45) 
over Q), the mode analysis is particulary simple since the re- 
sponse matrix C(ui) turns out to be already diagonal when 
the potential-density pairs (22, 21) are applied. This simplifies 
the calculations considerably, but puts severe limitations on 
the structure of the perturbations. All density profiles are real 
functions of the radius, excluding spiral structure. 

The distribution function for the Kalnajs disks becomes 
infinite at its outer edges. This singularity does not cause fun- 
damental problems in the mode analysis, but one still can ask 
to what extent the particular behaviour of these disks is de- 
termined by this feature. In order to address this question, we 
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finite, fitted to (45) and calculated the mode spectrum. Even 
for a relatively low value of g, the results are even quantita- 
tively very similar to those of the distribution (45). Both the 
mode frequencies and their spatial structure are in very good 
correspondence. As a sideproduct, this provides a convincing 
check for the validity of our calculations. 
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Fig. 7. Mode frequencies for model II with g = 3 and m = 2 

Fig 7 shows the mode frequencies as a function of the ro- 
tational parameter fi. In this case, the situation is much more 
complicated than in the simple diagram shown in fig. 3. It is 
not obvious anymore that the disk becomes more unstable for 
increasing fi, since there exist instabilities over the whole range 
of fi. In the case m = 0, there even exist unstable modes for 
fi = 0. Moreover, a particular unstable mode can become sta- 
ble again when fi is increased. When compared to fig. 7, the 
only striking common feature are the very rapidly growing un- 
stable modes occurring at high values for fi, which are actually 
the only ones that are present in Model I. This seems to indi- 
cate that these instabilities are a kind of "common" modes, in 
principle occurring in every uniformly rotating disk, regardless 
the structure of the unperturbed distribution. It is interesting 
to note that these "common" modes all appear close to co- 
rotation, having pattern speeds around fin, while most of the 
others have higher pattern speeds. In all the models which we 
analysed, these "common" modes were always present, while 



not a increasing function of the binding energy. The "com- 
mon" modes all appear at approximately the same value for 
Q, around 1.8. 

The response matrix C becomes more diagonal-like as the 
number of components used to fit (45) increases. Correspond- 
ingly, the imaginary part of the density profiles of unstable 
modes becomes smaller and smaller (in the limit of the Kalnajs 
models, the imaginary part vanishes), making the profile more 
bar-like and more aligned to the perturbing potential. 

2.3. Two-stream models 

In plasma physics, it is well known that two merged compo- 
nents flowing with different mean velocities can cause insta- 
bilities (the so called two-stream instability). One might ex- 
pect that a similar effect can occur in a galaxy consisting of 
multiple components with different streaming velocities (see 
e.g. Araki, 1987). Apart from the single-streaming Model I, we 
studied the stability of a galaxy resulting from the superpo- 
sition of two models I with a = 1 with equal but opposite 
rotation speeds 

df D ,o = tj(E* n +E*_ n ). (46) 

These distributions give rise to exactly the same mass density 
as in fig. 2. 



m = o 
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Fig. 8. Growth rate diagramme for Model I and two-stream Model 
I with a = 1 

When compared to the single-stream models, all modes are 
more stable except the lop-sided modes with m = 1 (the m = 
modes are not affected), which become less stable. This is in 
agreement with other investigations, both linear calculations 
(Araki, 1987) and numerical simulations (Zang & Hohl, 1978). 
This m = 1 unstable modes occur at Re(ui) = 0. Fig. 8 shows 
the stability diagramme for the m = 1 instabilities in the two- 
stream case, superimposed on that of the single component 
model. 

In the case of (46), the two counterrotating components are 
perfectly balanced and there is no net rotation. If one of the 
streams is weaker than the other, the m = 1 mode has some 
remaining pattern speed and shows a one-armed spiral shape, 
as shown in fig. 9. 

3. Conclusions 

We have explored the stability of stellar disks embedded in 
a quadratic potential. The main goal is to investigate the in- 
fluence of the distribution function on the stability. This has 




Fig. 9. Mass density of the lop-sided m = 1 instability in a partially 
counterrotating model I with a = and Q/Qn = 0.85. The net 
streaming velocity goes counterclockwise, i.e. the spiral structure is 
trailing. 



been achieved by studying the mode spectrum as a function 
of important parameters, such as slope of the binding energy 
dependence, velocity dispersion, mean rotation and disk/halo 
fraction. 

We have demonstrated that uniformly rotating disks are by 
no means the boring structures which they are often taken to 
be: they display a very varied stability behaviour, ranging from 
totally stable to violently unstable disks. Primarily responsible 
for this is the underlying distribution function. 

In our analysis, all isotropic non-rotating disks were sta- 
ble, if they have a distribution function which is increasing 
with increasing binding energy. The Kalnajs disks prove that 
the opposite is certainly not true. Although this does not prove 
any general behaviour, it might be that, at least for uniformly 
rotating disks, there is some law similar to Antonov's for the 
spherical isotropic case (Fridman & Polyachenko, 1984). As a 
general behaviour, non-inverted distributions become less sta- 
ble as they have more rotation velocity and less random ve- 
locity. They are stabilized again by increasing the amount of 
inert halo. 

An other remarkable property is the occurrence of spiral 
mode instabilities in most of the models. Usually, the forma- 
tion of spiral modes is associated with the presence of differ- 
ential rotation. It may well be that these modes can be seen 
as so-called "edge modes" (Toomre, 1981) or "groove modes" 
(Sellwood & Kahn, 1991), where amplification is triggered by 
the presence of a relatively steep gradient in mass density in 
the outer regions of the disk. As explained in section 2.1.2, the 
absence of spiral-like modes in the Kalnajs disks seems to be 
a mathematical coincidence rather than a general property of 
uniformly rotating disks. There are big qualitative and quanti- 
tative differences in mode frequency diagrammes between dif- 
ferent models, but they all tend to share a set of very rapidly 
growing modes if the disks are sufficiently cold. These com- 
mon instabilities seem to constitute the set of modes which 
are largely independent of the structure of the unperturbed 
distribution. 

The two-streaming version of model I is well suited to ex- 
amine the consequenses of the presence of two merged counter- 
rotating components on the stability. Of course, the stability 
analysis shows fundamental differences compared to the single- 



the pattern speeds. We confirm that the m = 2 and higher or- 
der modes are all considerably stabilized by the presence of the 
second stream. The m = is not affected, but the m = 1 mode 
is less stable than in the single model. 
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